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ABSTRACT 

We construct two-fluid equilibrium configurations for neutron stars with magnetic 
fields, using a self-consistent and nonlinear numerical approach. The two-fluid ap- 
proach — likely to be valid for large regions of all but the youngest NSs — provides 
us with a straightforward way to introduce stratification and allows for more realistic 
models than the ubiquitous barotropic assumption. In all our models the neutrons are 
modelled as a superfluid, whilst for the protons we consider two cases: one where they 
are a normal fluid and another where they form a type-II superconductor. We consider 
a variety of field configurations in the normal-proton case and purely toroidal fields 
in the superconducting case. We flnd that stratification allows for a stronger toroidal 
component in mixed-field configurations, though the poloidal component remains the 
largest in all our models. We provide quantitative results for magnetic ellipticities of 
NSs, both in the normal- and superconducting-proton cases. 
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1 INTRODUCTION 

It has long been acknowledged that neutron stars (NSs) have strong magnetic fields, which may play a variety of roles in the 
physics of these objects. Magnetic fields can change the stability of a NS, affect its evolution and cooling prope rties, induce 
precession and provide the restoring force for a class of oscillation modes ( Harding fc^^aj 200^ Mereghetti 2008h. They ma y 
also produce crust-core coupling and affect the post-glitch spin behaviour of NSs ( Easson 19791 : Dib. Kaspi fc Gavriil 20081 ) . 
All of these effects are likely to be sensitive to the details of the magnetic-field configuration in the stellar interior, whereas 
we are only able to observe the exterior, dominantly dipolar, field. For this reason considerable effort has gone into modelling 
the interior magnetic fields of neutron stars; we give an overview of this in section 2. 

Despite considerable progress, the vast majority of studies on magnetised neutron stars treat them as barotropic single- 
fluid bodies. Whilst this is a sensible way to begin producing simplified models, more realistic studies need to confront 
some of the detailed interior physics of NSs. For example, over forty years si nce the first predictions that neutron star 
matter would contain superfluid neutrons and superconducting protons (see, e.g., Bavm. Pethick fc Pines ( 19691 )). there have 
been no attempts to construct multifluid magnetic equilibria along these lines. The issue is highly topical, too — recent 
observations of the cooli ng of the young NS in Cassiopeia A are well explained by a model with superfluid neutrons and 
superconducting protons ( Shternin et al.l201ll : Page et aPboilh . In the case of no magnetic field there h as been more progress, 
includ i ng studies of the equi l ibria and oscillation modes of superfluid NSs (e.g . Prix fc RieutordI (2002); Prix. Novak fc Comer 
lYoshida fc Eriguchil (|2004l ) : iPassamonti. Haskell fc Anderssonl ((20091)). The closest related work on magnetised stars 
is probably the study of smgle-flmd barotropic NSs with superconductivity, bv lAkgiin fc WassermanI (|2008l ). 

A simpler problem than a NS with superfluidity /superconductiv ity is a stratified (i.e. non-baro tropic) , magnetised star 
composed of normal fluid. Stratification is expected to exist in NSs (jReisenegger fc Goldreichlll992l ). where it can allow for 
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a wider range of equilibrium configurations as well as providing a stabilising effect o n the magnet ic field ( Reiseneggei 20091 ) . 
Again, the limitations of barotropic stellar models have been known for decades (Mestel 19561 ). but there have been few 
attempts to construct stratified magnetised equilibria. Most of the work in this direction has been that of Braithwaite and 
collaborators, who perform nonlinear magnet ohydrodynamics (MEID) evolut i ons of a stratified star and show that the system 



appea rs to relax to an equilibrium over time (jBraithwaite fc Nordlundlliooel : lBraithwaitell2009l ). In addition, iMastrano at al. 



very recently studied single-fluid non-barotropic models of NSs, albeit using pre-specified magnetic field configurations. 

In this paper we attempt to plug some of the gaps in the literature, presenting the first results for magnetised NS equilibria 
modelled as a two-fluid system. We treat the neutrons as being a superfluid, which mainly interacts with the proton fluid 
through the star's gravitational potential. By choosing different degrees of compressibility for the two fluid species, we can 
introduce stratiflcation. We begin by considering the case where protons are a magnetised normal fluid and find equilibria with 
purely poloidal, purely toroidal and mixed-field configurations. After this, we treat the 'full' problem where the neutron star 
is composed of superfluid neutrons and type-II superconducting protons for the first time, specialising to the case of a purely 
toroidal magnetic field. This paper takes a nonlinear approach to the problem, but we have also worked on a separate st udy 
using perturbation theory; the results of this are presented in a parallel paper. iGlampedakis. Andersson fc Lander! (|2011f ). 

This paper is laid out as follows. In section 2 we present the equilibrium equations for our stellar model, independent of 
the details of the magnetic field. We then specialise to the case where the protons are a normal fluid governed by standard 
MHD. In section 3 we consider the case where the protons are instead a type-II superconductor, presenting simplified equations 
which nonetheless are likely to be accurate enough for an initial study of equilibria. In the case of purely toroidal fields we 
derive equations in a form suitable for numerical integration. In section 4 we give details of our numerical procedure, whilst 
in section 5 we show how to convert code output (in terms of dimensionless code variables) into useful physical quantities. 
Section 6 contains our results for two-fiuid equilibria, exploring the effect of stratification and type-II superconductivity and 
section 7 is a discussion. 



2 TWO-FLUID EQUILIBRIUM EQUATIONS AND NORMAL MHD 



2.1 Equilibrium equations for two-fluid magnetic stars 

We begin by modelling a neutron star as axisymmetric, with the magnetic symmetry axis being the same as the rotation axis. 
For this system it is natural to use cylindrical polar coordinates {zo, 4>, z), where the z-axis is aligned with the symmetry axis 
of the star. Although relativistic effects will certainly play a role in neutron star physics, we will neglect these for this study, 
and work in Newtonian gravity. We also ignore the elastic crust and any non-hadronic matter that may form an inner core, 
and assume that the neutron star is a multifluid system, with a neutral fluid of neutrons and a charged fluid of protons and 
electrons. 

Our multifluid model should be applicable for the bulk of a neutron star's interior during much, but not all, of its life. A 
neutron star is born hot, with its charged and neutral components strongly coupled through various dissipative mechanisms; 
at this stage it behaves as a single, non-barotropic, fiuid. When the star cools below a critical temperature Tc however, the 
core neutrons begin to condense into a superfluid stats0. The recent observations of the cooling of the Cassio peia A NS 
indicate that Tc ^ (5 — 9) x 10* K; a NS core should start to drop below this temperature ~ 100 years after birth ( Page et al. 
20 111 : IShternin et aDboill ) . For this reason, even young pulsars are likely to hav e substantial muhiflui d regions — i.e., regions 



where the neutrons are superfluid. Models of magnetar temperature profiles ( Kaminker et al. 20061 ) indicate they are cool 
enough to have multifluid regions too. For this preliminary study into stratified NSs we model the entire star as multifluid; 
this is obviously a simpliflcation, and more detailed modelling should account for the layered structure (including single-fluid 
regions) expected in real NSs. 

Now, since electrons have a much smaller mass than either of the other components, we take the standard approach of 
neglecting their inertia. Our final model is then a two-fiuid system, with a proton fiuid (whose variables are denoted with 
a subscript roman index p) and a neutron fluid (with index n). Neutrons and protons are assumed to be of equal mass 
rrin = rrip = m. Next, we denote the particle chemical potentials by /^n and /ip, then define /in = Hn/m and fi^ = {fip -\- ij,c)/m\ 
note that we use fi^ to denote the chemical potential of the whole charged component, as it also contains a contribution from 
the electrons /io. 

In normal MHD, only the protons are coupled to the m agnetic field, but for a superfiuid /superconducting system t here 
is generally a magnetic force acting on the neutrons as well ( Mendelll 1991 ; Glampedakis. Andersson fc SamuelssonI 201lh . In 
this work however, we only consider the case without entrainment, for which the neutr on superfluid is decoupled from the 
magnetic field in the superconducting case too. This issue is discussed in more detail in IGlampedakis. Andersson fc Lander 



^ Neutron superfluidity in the inner crust is likely to occur at a higher temperature — above 10 K — due to its singlet pairing type 
(as opposed to the triplet type in the core), making this part of the star a multifluid system too. 
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( 2011f ). Given these assumptions, we may now state the separate Euler equations which govern the two fluids: 

v('/i„+4-4SUo. (1) 



where Fmag denotes the magnetic force; the form of this force will be the only difi'erence between the normal-MHD and 
superconducting cases to our order of working. As usual, $ denotes the gravitational potential, p is density and fl denotes 
rotation rate (with indices for individual fluid species). We will only consider the case where there is no entrainment and the 
two fluids corotate, i.e. fin = ilp = fi. Since we are only studying stationary configurations, the equations do not contain 
mutual friction terms. 

As in the single-fiuid case we have Poisson's equation for the gravitational potential, from which we see that the behaviour 
of the two fluids is linked despite the neutrons being a superfluid: 

V'$ = 47rGp = 47rG(pn + Pp). (3) 

Let us multiply equation ([T| by pn, equation ([2]) by pp, and add them: 

P,iV/in + PpV/ip + pV$ - pV ^ ^ = Fmag. (4) 



2 

Now since pnVjln + ppVjlp = VP (the gradient of the total fluid pressure), this allows us to recover the familiar 'single-fluid' 
Euler equation: 

^+v$-vf^U^:i^, (5) 
p \ ^ J p 

except that in this two-fluid case one will typically have stratiflcation, and hence can no longer replace the pressure term 
with an enthalpy gradient. Instead, it is more useful to work with the proton-fluid Euler and a 'difference-Euler', given by ^ 
minus fTJ: 

V(/ip-/In) = ^^. (6) 

Taking the curl of this equation, we see that (as in the single-fluid case) there exists a scalar M such that 

5^ = VM. (7) 

pp 

The magnetic force — whether the normal MHD 'Lorentz force' or the flux-tube tension force of a type-II superconductor 
(see section 3) — therefore depends on a scalar function M and the proton-fluid density pp, as opposed to the single-fluid case 
where the dependence is on the total mass density: Fmag/p = VM. In deriving the equations that govern the behaviour of 
the magnetic fleld in the two-fluid case, the only differences from the single-fluid case will be in this density factor; the flnal 
magnetic equations will be the same but with p replaced by pp. 

For our numerical scheme we need to work with the equilibrium equations in integral form; in this section we discuss 
those equations which hold for all magnetic-fleld conflgurations. Firstly, we always have the same form of Poisson's equation 
for the gravitational potential: 

^_^r Pn(O + Pp(r0 

./ r - r' 



Pp+$-^ = Af + Cp (9) 



We also have flrst-integral forms of the proton-fluid Euler equation ([2]) 
and the difference-Euler ([6]): 

Pp - Pn = M + Cd, (10) 

where Cp and Cd are the integration constants for the proton and difference-Euler equations, respectively. Finally, for the 
equation of state we choose an energy functional £ given by 

£• = fcnp?' fcpp^^ (11) 

a two-fluid analogue of a polytropic model ( Prix. Comer fc Andersson 20021 : Passamonti. Haskell fc Andersson 20091 ). This 
could also be written in terms of the polytropic indices A'x = l/(7x — 1), x = {n, p}. Note that can be generalised to 
include cross-terms, corresponding to symmetry energy and entrainment. The chemical potentials are then deflned from the 
energy functional by 

-^£1 • 
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where the index x represents one particle species (neutrons or protons) and the y index represents the other. The other 
equations of the system depend on the details of the magnetic-field configuration; we discuss these next. 



2.2 Normal MHD 



The simplest case is when the charged component of the stellar matter is governed by normal MHD. The magnetic force Fmag 
is then the familiar Lorentz force, given by 

1 



£ = j X B = —(V X B) X B, 

47r 



(13) 



where B is the magnetic field and j the current. Most studies of NS equilibria have used the normal-MHD equations and have 
additionally adopted a single-fluid model, whe re the charged component is t he entire fluid — this an s atz has been used b y 
innumerable authors. R ece nt examples include Tomimura fc Eriguchi (2005); Kiuchi &c Yoshida (2008); Haskell et al. ( 20081 ) : 
Lander fc Jones (200£) and Ciolfi. Ferrari fc Gualtieril ( 201C ). but there are decades of earlier work on single- fluid models — 



see, e.g.. 



Ferrarol (|l954[ ) and lRoxburghl (jlQed ) 



An obvious way to progress is to work with a two-fluid model, but still in normal MHD: the protons are a nor- 
mal fluid and the neutrons su perfluid. This is not just a stepping stone towards the superfluid/superconducting model of 
Bavm. Pethick fc Pine^ (|l969h — it already allows us to study the role of stratification in magnetised stars. Furthermore, 



the interior magnetic field strength in magnetars may well exceed the second critical field Hc2 ~ 10^^ G at which supercon- 
ductivity is destroyed; in this case magne tar matter could indeed be predominantly normal protons and superfiuid neutrons 
I Glampedakis. Andersson fc LandenuOllf ). 

Since only the proton fiuid is magnetised, the derivations for all equations with magnetic terms proceed in the same 
manner as the single-fluid case, with factors of p replaced by pp; otherwise they are identical. For this reaso n we do not give 
the de tails here, but the essentially identical single-fluid derivation may be found in various papers (see, e.g.. lLander fc Jones 
The non-magnetic equations of the system are those discussed in the previous subsection. 

In the mixed-fleld (or purely poloidal field) case we have a Poisson equation for the ^-component of the magnetic vector 
potential too, allowing for the magnetic field to extend outside the star: 

lfiu)f'{u)+47T^M'{u)pp{r 



A4,(r)sin(j} = — 



sin cj> df , 



(14) 



where u is a streamfunction; see the next section. In this equation / and M are functions of u, tildes denote dummy variables 
and primes denote derivatives. In addition, u and are related through u — u^A^, so the integral in equation (|14p may be 
rewritten in terms of A^. As mentioned earlier, this integral equation for A^ is identical to the single-fiuid version except with 
a p term in the integrand replaced by pp. 

The magnetic functions M{u) and /(«) appear to be arbitrary, but in practice we h ave found only o ne acceptable 
functional form for each; this was also the case in our earlier study of single-fluid models ( Lander fc Jonej 2003 ) . In the 
mixed-field case, we use the following form of M: 

Mmix = KU — HVjA^ (15) 

where is a constant governing the strength of the Lorentz force. We flnd that choosing M{u) as a lower or higher power of 
It results in the code, through equation (|14p . iterating to an unmagnetised equilibrium solution — and hence only the above 
choice of Mmix is acceptable in our work. 

The other magnetic function / in the integral equation for relates to the toroidal component of the fleld and is defined 
to be non-zero only within the largest interior closed field line, so there is no exterior toroidal field (and hence no exterior 
current): 



f{u) = a{u - Umax) ' H {u - Umax), 



(16) 



where Umax is the maximum value attained by u within the star and H is the Heaviside step function. The strength of the 
toroidal component is adjusted by varying the constant a. Other choices of the exponent than 1.1 are possible; we choose 
this value as it gives a relatively strong t oroidal-field componen t. Choosing a value of unity or less for the exponent gives 
undesirable behaviour in f'{u), however ( Lander fc Jonej 200^). Finally, whilst adjusting the exponent in (|16p allows for 
configurations with very modest differences, we do not know of any genuinely different functional form of / which is still 
dependent on u and restricts the toroidal component to the stellar interior. 

In the toroidal-field case there is no separate integral equation for the magnetic field. M has a different functional form in 
the pure-toroidal field case Mtor from that in the mixed/pure-poloidal field case Mmix- As with other quantities, the two-fiuid 
form of Mtor may be found from the single fiuid version by replacing p with pp-. 



Mtor 



1 

47r 



h{C) dh 



dC, 



(17) 
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where h is a function of Pp'cu^, related to the magnetic field by 



(18) 



We choose h{ppzu^) = Xp^m^ , where A is a constant governing tlie strength of the toroidal field. Once again, all other functional 
forms appear to result in either divergent magnetic quantities or our numerical scheme iterating to an unmagnetised solution 
( Lander fc Jonedboogl V 

Note that the difference-Euler ((6| for an unmagnetised star shows that the two fluids are in chemical equilibrium. In 
the magnetic ca se there is an extra Fmag term — this c ould be interpreted as a statement that the star is out of chemi- 
cal equilibrium (jGlampedakis. Andersson fc Landeill201ll '). Alternatively, one could adjust the notion of what is meant by 
' chemical equilibrium'. Since the energy func tional £ gains an extra contribution from the presence of a magnetic field 
I Glampedakis. Andersson fc Samuelssonll201ll ). the chemical potentials — defined through (|12p — also change accordingly. 
One could then see equation (|6]l as quantifying the difference between magnetised and unmagnetised chemical equilibria. 



3 MAGNETIC FORCES IN TYPE-II SUPERCONDUCTIVITY 



In a typical N S the protons are likely to form a type-II superconductor, where the field penetrates the star through 
thin fiuxtubes ( Bavm. Pethick fc PinesI 1969h . This results in a magnetic force de pendent on the averaged magnetic field 



B but also on the first critical field Hci- T he result is a fluxtube tension force (jEasson fc Pethickl 119771 : iMendelll Il991 



Glampedakis. Andersson fc Samuelssonll201ll '). quite different from the Lorentz force (|TS 



The aim of this section is not to give a complete description of the equations governing a star with type-II supercon- 
ducting protons , but i nstead to explore how many of the steps from the normal-fluid Grad-Shafranov derivation (see, e.g., 
Lander fc JonesI |2009h l also hold in this case. We derive a general 'interim' expression for the magnetic held, analogous 



to the normal-MHD case, at which point one has to specialise to either purely toroidal flelds or mixed flelds (with purely 
poloidal flelds as a particular limiting case). A key result for mixed fields from the normal-MHD case no longer holds in 
the superconducting case, and we defer a detailed study of this to a later paper. Instead we concentrate on purely toroidal 
magnetic fields, where the derivation is more straightforward. 



3.1 General form 



In iGlampedakis. Andersson fc SamuelssonI (|2011f ). a full derivation of the equations governing a type-II superconducting 
neutron star is presented. For our purposes, however, only a reduced and simplified set are needed. We consider the case 
where there is no magnetic force on the neutrons, and neglect small terms like the 'London field'. Given this, our equations 
reduce to those of section 2.1, by setting Fmag to be the superconducting magnetic force Smag- Before discussing this force, 
we have one useful result from normal MHD which still holds in this case. Usin g the fact that V ■ B = together with the 
assumption of axisymmetry allows us to write B in terms of a streamfunction u ( Lander fc Jonej 20091 '): 



B 



B 



tor = — Vit X ej, + B^es- 

VJ 



(19) 



As in normal MHD then, we have B- = 0. At this point we turn to equation (153) from lGlampedakis. Andersson fc Samuelsson 
( 20111 ) for the explicit form of Smag^ 



^mag = ^[(B-V)(/f,lB) 



(20) 



where B is a magnetic unit vector and B is the (position- dependent) mag nitude of the field, so that B = SB. From the 
calculation in appendix A2 of [Glamped akis. Andersson fc Sa muelssonI (|2011 ). the first critical field Hd ~ hcPp/e* to a good 
approximation, where e* is entrainment and he a constant parameter. Since we are assuming that there is no neutron-fluid 
magnetic force, and hence no entrainment, we have = 1. The magnetic force then becomes 

Smag = ^ [(B • V)(ppB) - V(Bpp)] . (21) 
Now by using a vector identity to rewrite the flrst term on the RHS and rearranging the result we arrive at: 

„X ^22) 

V X B; 

(23) 



Next we want to rewrite the bracketed term on the RHS of equation 
as in the single-fluid case it may be shown that 

J = — V(ci7Brf,) X e^ -f uerf,. 



PpVB + B X ( ppV X B + Vpp X Bj . 

Let us start by deflning a 'unit current' j 
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We then decompose both of the bracketed terms from (|22|) into poloidal and toroidal components. The resuhant expression 
may then be rearranged to show that 

PpV X B + Vpp X B = —V{ppwB^) X + ( ppj^ - YfP-Y^\ (24) 



The magnetic force becomes 



r-5mag = Pv^B + — B X ( V(/9pII7B0) X + ( PpJ> - ) B X G^. (25) 



he VJ \ / \ VJ±> 

Again, it is useful to rearrange some terms in this expression. The second term on the RHS of (|25p may be rewritten as: 

— B X (v{ppUjB^) X e^) = -\vu X V{ppmB^) + ^V{ppZuB^), (26) 
using standard vector identities. To simplify the third term on the RHS of (|25p we use 



B X 60 = Bpol X = . (27) 

tx7 



Now, putting these two results into the expression for the magnetic force ([25} we find that: 

- I^Smag = PpVB + -\vu X VippmB^,) + ^VippZuB^,) + [ Z^P^ _ PpJ^ ) Z!^. (28) 

tic ^ ^ \ ttjij / UJ 

By axisymmetry, the toroidal component of S^ag must be zero. Now, all of the terms in the above expression are gradients 
of scalars (and hence poloidal), except the cross product. This term is purely toroidal and hence must be zero: 

AtVm X V{ppZuB4,) = 0. (29) 

Hence we can remove this term from the magnetic force to arrive at a somewhat simplified result, in terms of the gradients 
of various scalar functions: 

47r^ _ 47r ^ B^^^ h \ , (^Pv ' . ^ ^ 



l-5mag = -7-ppVM = ppVB + ^VippmB^,) + PpJ0 , (30) 

he he taj \ run J zu 

where M is defined by equation as before. For compariso n, the equivalent form for the magnetic force in the normal- MHD 



where M is deimed by equation (IVj), as beiore. i^or compariso n, t 
derivation (see, e.g., the appendix of iLander fc Jones (|2009l )') is 



AnpVM = ^V{vjB^) - 4nj^— (31) 
zu zu 



where = j^[V x B]^. We see that in the superconducting case there are extra terms with VB and Vpp factors^ 

The above derivation is similar in approach to that for the Grad-Shafranov equation of barotropic normal MHD ( Lander fc Jonej 



20091 ). up until the result H30p . One remaining step from the normal-MHD derivation cannot, however, be generalised straight 



forwardly for superconducting matter: we no longer have B • VM = and so in general M 7^ M{u). This prevents us from 
continuing the derivation for poloidal/mixed fields in the normal-MHD manner. Rather than discussing this issue in more 
detail now, we stop at the interim general result pO|l and only consider the simplest specific case, where the magnetic field is 
purely toroidal. We intend to return to the cases of purely poloidal and mixed poloidal-toroidal fields in future work. 

3.2 Purely toroidal fields 

For purely toroidal fields, let us return to equation pop . In this case we have Vu = and the expression for the magnetic 
force reduces to 

Smag = Pp'^M = -!hLlv{zuppB^), (32) 
47r zu 



where we have also used B^ — 1. Again, we recall the normal-MHD result ( Lander fc Jones 20091 ) for comparison 



C = pVM = -^^V{zuB4,). (33) 
Now, dividing equation (|32p by pp and taking the curl of it yields 

= V ( — !— ) x \/{zuppB4,) = l—V{zupp) X V{zuppB^). (34) 

\ZXjppJ uj^p^ 

Wherever the magnetic field is non-zero zupp 7^ too, since this is the geometry of a toroidal field, so we have 

= V(c;7pp) X V{vjppB^) (35) 
and hence the two arguments of the gradient operators are related by some function r]: 

Vizupp) = TXjppB^. (36) 

Putting this into (|32p and defining = zupp we get 

„ , , he 1 „ , he 1 Ari „ ^ ,„„. 
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where the corresponding magnetic field is 



B = B^e^ = ^e^. (38) 



The first integral of (|37[) gives our final result, which may be used in the code: 

Not all functional forms of 77(^) are acceptable choices — for example, taking 77 to be a linear function of ^ leads to a magnetic 
force which diverges on the polar axis. In this paper we will work with two other choices of rj. The first (which we will refer 
to as 'C^-superconductivity' for brevity) is 

V = '^oC' (40) 
with rjo a constant, which may be varied to adjust the magnetic field strength. With this we have 

M = — ^'^''^ ropp and B4, = rjo-oopp. (41) 

This form of is the same as that in our normal-MHD toroida l-field case. Note that ("^-superconductivity is the two-fluid 
equivalent of the case considered by Akgiin fc Wasserman ( 20081 ) . The second functional form we will work with (hereafter 
'C^-superconductivity') is 

V = %C'- (42) 
The magnetic force scalar and magnetic field in this case are given by 

M = -^^^vj^pI and B^^r,ovj^pl. (43) 

At this stage, our motivation for choosing these two functional forms is their simplicity. Later, we will find that the resulting 
equilibria in the two cases are very similar, suggesting that our results may be quite generic. 



4 NUMERICS 



4.1 Overview of code 

The code we use iteratively solves the equilibrium equations for our neutron star model, and being non-linear is not restricted 
to the perturbative regime of slow rotation and weak magnetic fields. Our numerical scheme is b ased on the Hachisu sel f- 
consistent field (SCF) method ( Hachisu 19861 ). a more robust extension of an earlier SCF method by Ostriker fc Mark ( 19681 ). 
A version of the Hachisu SCF method for magnetised stars was presented by Tomimura fc Eriguchil ( 2005 ). Since our scheme 



is a fairly straightforward extension of these p revious ones, w e co ntent ourselves with a summary here and focus on the 
differences between them; we refer the reader to Hachisu ( 19861 ) and Tomimura fc Eriguchi ( 2005 ) for more details. 

To find equilibrium models, the user must specify a number of stellar parameters at the outset. The major ones are 
related to: 

1. the EOS — through the polytropic indices A^n and Np (related to the compressibility of each fiuid); 

2. the shape of the star (related to the rotation rate) — specified through the ratio of polar to equatorial radii rpoie/r^q, and 
the ratio of the neutron-fluid equatorial surface to the proton-fluid equatorial surface r^g/r^g; 

3. the magnetic field — through k and a for normal-MHD mixed fields (or poloidal fields) , A for normal-MHD toroidal fields, 
and 770 for toroidal fields in a superconductor. See sections 2.2 and 3.2 for more details. 

We nondimensionalise all quantities within the code using the requisite combination of r^q, G and pmax and use a hat 
to denote these dimensionless quantities. As a consequence we have req = Pmax = 1. When calculating integrals we first 
decompose the integrands into radial pieces and Legendre polynomials (note that there is no azimuthal piece, since we work 
in axisymmetry). The contributions over all relevant grid points are then summed up appropriately. We include angular 
contributions up to degree / = 32 in the Legend re polynomials, thus allowing for magnetic configurations of high multipolar 
structure, in contrast with the dipolar models of iGlampedakis. Andersson fc Landed (|201ll ). 



4.2 Finding integration constants and the rotation rate 

To find the stellar rotation rate and integration constants of the Euler equations we have to evaluate these equations at suitable 
boundaries. Typically the two fiuid surfaces do not coincide, and there are different procedures to calculate these constants 
depending on whether the proton or neutron fluid is outermost. However, we know a posteriori that all our configurations 
have the proton fiuid either outermost or coincident with the neutron surface, so we only consider that case here. 

For this, we work with the proton-fiuid Euler and the difference Euler and aim to find Cp, d and Q.; recall that we assume 
that both fiuids have the same rotation rate. Sin = f^p = In this section we define each fiuid surface by the vanishing of 
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its corresponding chemical potential /ix = 0; this approach was also used bv lYoshida fc Eriguchil (|2004 ) and avoids numerical 
difficulties associated with dealing with the true fl uid surfaces — those given b y px = 0. In general the two deflnitions do not 
agree, due to coupling between the different fluids ( Prix. Novak fc Comer 20051 '). and in practice one would expect an interior 
phase transition rather than a sharp fluid surface. Despite these issues, working with the /ix = surfaces is sufficient for our 
simplifled model of a two-fluid NS. In section 5.2 we describe how to check the surfaces of the redimensionalised, physical, 
star. 

This subsection requires the user-specified stellar axis ratio rpoie/req (which is the same as the proton-fluid axis ratio 
^poie/^eq)- In our dimensiouless units the axis ratio is equal to the dimensionless polar radius of the proton fluid r^^j^. The 
neutron-surface equatorial radius f"^ — r^^q/^eq is also used here. We begin by evaluating ^ at the equatorial stellar surface 
feo = 1: 



Ap(l) = = C7p + M(l) + — -$(1). 
and at the polar stellar surface rpoie{~ r^^i^): 

Mr pole) = = Cp -f M(fp) - $(fp). 
We may combine the above expressions to find Q, and Cpi 

= 2(M(fp„,e) - Af(l) - ^fpoie) + -I'(l)) 



(44) 
(45) 

(46) 
(47) 



Note that we have M oc pp for a toroidal field (in both normal and superconducting matter) and so is zero at the stellar 
surface; in this case the expressions for Q, and Cp have no explicit dependence on the magnetic force. Now evaluating the 
difference Euler at the neutron-fluid surface, rather than the stellar surface, we determine Cd- 

Cd = Ap (r^q) - M {r:^) . (48) 



4.3 Iterative step 

We recall that the chemical potentials are defined from the energy functional £ by 

Ax = 1^ . (49) 
For our particular functional the resulting expressions for the chemical potentials may be rearranged to give 

Pn=(^) ,Pp=f^) . (50) 

Evaluating these at the centre of the star we obtain relations between the maximum values of the chemical potential and 
density for each fiuid. In dimensionless form these may be rearranged to show that 

fcp7p = f^r^'^AO)-'^"'' , fcn7n = Ar''(l - X-p(O))-^/^", (51) 

where x-p(O) is the central proton fraction, using the fact that pp — Pp/p = Xp and similarly pn = (1 — Xp). Comparing 
these relations with those in (|50l) we obtain expressions which will allow us to iterate for new density distributions within our 
numerical scheme: 

/ - \ JVp / - \ iV„ 

Pp--p(0)(^) , An = (l-.p(0))(^j . (52) 

4.4 Iterative procedure 

We now have all equations needed for the iterative procedure that generates our equilibria. As discussed at t he start of this 
section, the procedure is an extension of the Hachisu SCF method (jHachisulllQsi : iTomimura fc Eriguchill2005l ). To begin the 



iterative process (step 0), we set the two fluid densities and the magnetic potential to be constant within the star; however, 
we have found that the final equilibria are independent of the values chosen at this step. Note that for a toroidal field — in 
normal or superconducting matter — there is no separate iteration for the magnetic field and hence no step 2; instead the 
behaviour of the field is directly linked to pp. The iterative steps are as follows: 

0. For the initial iteration, start with a guess that pn, Pp and are all constant; 

1. Calculate the gravitational potential <1> from the pn and pp distributions and Poisson's equation ((Sj; 

2. Calculate the new magnetic potential component from its value at the previous iteration A^''', using the magnetic 
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Poisson equation (|14p with A^'^ and pp in the integrand; 

3. Evaluate the proton-fluid Euler ((9)1 at the equatorial and polar surface, using boundary conditions on flp; this gives two 
equations which fix and then Cp; 

4. We are now able to use the proton-fiuid Euler to find fip throughout the star; 

5. Evaluate the difference-Euler (|10p at the equatorial neutron-fluid surface to find Cd', 

6. Now use the difference-Euler to find /in throughout the star; 

7. Calculate the new density distributions from the chemical potentials, using the expressions in (|52|) : 

8. Return to step 1 using the new Pn,Pp and A^; repeat procedure until satisfactory convergence is achieved (i.e. until the 
difference between quantities at consecutive iterative steps is less than some small tolerance). 



5 STELLAR PARAMETERS AND SEQUENCES 
5.1 Redimensionalising 

When looking at sequences of configurations, to study the effect of magnetic fields or rotation, we need to ensure we are 
always working with the same physical star — the dimensionless quantities produced by the code are not sufficient. Two stars 
are physically the same if they have the same mass and equation of state: M,^n,^p,kn and kp must all be equal for both. 
Since we specify the (dimensionless) indices 7n and 7p at the outset, we only need to fix the physical values of A^, fcn and hp. 

Fixing the stellar mass is straightforward — we choose Ai = IAMq (where A^0 is one solar mass) — but more care 
is needed with the polytropic constants. We choose to fix their physical values in the following way: take a spherical star 
in hydrostatic equilibrium (nonrotating and with no magnetic field), with some fixed 7n and 7p and coincident neutron and 
proton surfaces. Now find the polytropic constants fen, fep that give it a radius of 10 km. More specifically, the dimensionless 
mass and polytropic constants are related to their physical counterparts through: 

M f fcx 



M = 3 , fe"x = — 2 . (53) 

Pmaxrg, Grig 

Having specifying a spherical configuration with some particular 7x, the code may be used to calculate M and fex. The 
requirement that A4 — 1.4A^0 and Veq = 10 km for the spherical configuration then fixes the physical value of fex for the 
stellar sequence; from this we may find pmax and r^q for any star in the sequence, and hence redimensionalise any other 
quantities (e.g. rotation rates and field strengths). Any rotating, magnetised configuration redimensionalised in this manner 
will be the same physical star: if you stopped it from rotating and removed its magnetic field, it would return to the same 10 
km spherical body. 



5.2 Adjusting the position of the fluid surfaces 

Our iterative procedure does not automatically adjust the location of the neutron surface, even though it will change with 
respect to the proton surface when the star is rotating or magnetised. A rotating or magnetised star in which the two fluids 
coincide at the equator would not return to our canonical unmagnetised nonrotating spherical star with coincident fluid 
surfaces. This means that the physical values of the polytropic constants fex would be different in the two cases. Instead we 
have to manually adjust the location of the neutron-ffuid equatorial surface so that the values of fcx do agree with those of 
our canonical physical star. In practice the difference between the fluid surfaces is only significant at high rotation rates or 
magnetic fleld strengths, and for the purposes of this study a simple root-finder algorithm can be employed to run the code 
a few times, adjusting the neutron-surface until satisfactory results are achieved. 



5.3 Virial test 

The scalar virial theorem, which is derived from the (sum) Euler equation, states that a certain combination of energies 
balances the acceleration of the system. For a single-fluid polytrope: 

2r + £:n.ag + lV + 3(7-l)f/= i0 (54) 

where T, fmag, and U are, respectively, the kinetic, magnetic, gravitational and internal energies of the system, and / is 
the moment of inertia. For our 'separable' equation of state, where each fluid obeys its own polytropic relation, this relation 
may be generalised simply: 

2T + £:„,ag + 3(7„ - 1) [/„ -f 3(7p - 1) [/„ = i § . (55) 

2 at'' 

Since we seek equilibrium configurations, we can test the accuracy of our code by how close the residual acceleration of the 
system (i.e. the right-hand side of the previous equation) is to zero. We then form a dimensionless quantity by normalising 
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this acceleration through division by W: 



Virial test' . error = 1^^ + + + 3(7n - 1)^. + 3(> - l)^p| 

\W\ 



5.4 Ellipticity 

Our numerical scheme involves specifying some axis ratio rpoie/rEq, corresponding to the surface shape of the star. A more 
informative quantity however is the ellipticity, which measures the deviation of the whole mass distribution from sphericity. 
Starting from the mass quadrupole moment tensor 

Qjk = j pXjXk AV (57) 

where Xj , 3;^ are Cartesian coordinates, we define the ellipticity of the star through the components of the quadrupole moment 
at the equator Qxx and pole Qzz'- 



e = 




5.5 Magnetic-field quantities 

The conventional definition of the (normal-MHD) magnetic energy is as an integral out to infinite radius: 

1 



57r 

all space 



(59) 



which would have to be truncated on our finite numerical grid. One alternative but equivalent form o f £mag is as the integral of 
r ■ C; this may be seen from the relevant part of the virial theorem derivation (|Chandrasekhailll96ll ) . Since C is only non-zero 
over the star (by virtue of its pp factor) , we may write the magnetic energy as 



v-CAV. 



The magnetic energy in the toroidal component of the field is also an integral over the star: 



Bl AV, 



(60) 



(61) 



since = outside the star. Finally, it is clear that the poloidal-field contribution to the energy is simply £mag — £tor, so 
there is no need to use an integral out to infinity for this case either. 

When looking at mixed-field configurations, we want some measure of the relative strength of the toroidal and poloidal 
field components. We do this in two ways: one related to the global contributions of each component and one which looks 
at their maximum values. For the former, we take the percentage of the magnetic energy stored in the toroidal component 
ftor/£mag; this is relevant for understanding quantities like the ellipticity, which scale with . For the latter we compare the 
maximum value attained by each field component, using the ratio 

-"tor 



Dma 
-'^tor 



I Dmax 
+ -"pol 



(62) 



Unlike the more obvious choice of B^f^ /B^^^, the ratio (|62p varies from zero (no toroidal field) to unity (purely toroidal 
field), which is perhaps clearer. 

Alth ough a number of s tudies have used the energy ratio f^tor/fmag to assess the stability of a magnetic-field configuration 
(see, e.g., Braithwaitel ( 20091 ) ). thi s is a global quantity, w hereas the relevant instability is highly localised — around the neutral 
point for a purely poloidal field ( Markev fc Tavler 19731 ). To remove this in stability, wha t is more important is whether the 
toroidal component is locally comparable in strength with the poloidal one (| Wright! 1 197^ 1. Since the strength of the poloidal 
component in the closed-field line region is directly related to its maximum value, we argue that the ratio (|62[) may be a 
better indication of the stability of a configuration. 

Since we are interested in the effect of a magnetic field on the global properties of neutron stars, we work with a volume- 
averaged form of the magnetic field B rather than a surface value: 

1 f „2 STTfnrag ^gg^ 



B^ AV 



V J V ' 

all space 

where V is the star's volume. This volume-averaged value is approximately double the polar surface value Bpoie for a poloidal 
field (for a toroidal field the surface value is zero). 

For superconducting configurations we have to fix the value of the first critical field strength Hd- This is not a constant. 
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virial test 




grid points 

Figure 1. Testing the order of convergence for our code, by plotting the virial test result against the number of grid points. All results 
arc for a rotating magnetised NS model, stratified with Af,i = 1 and A^p = 2, with B = 1.2 X 10^^ G and Q = 650 Hz. The line drawn 
shows the expected result for second-order convergence; it agrees well with the actual points from different resolutions. 



Table 1. Comparing with the unmagnetised two-fiuid models presented in Prix and Rieutord. The values shown for the polytropic 
constants and masses arc our results, in our dimensionless units, with the discrepancy from those of Prix and Rieutord indicated in 
brackets. 





7n 


7p 




kp 


M 


Model I 


2.0 


2.0 


0.7074 (0.1%) 


6.366 (0.1%) 


1.273 (0.5%) 


Model II 


2.5 


2.1 


0.7060 (0.5%) 


9.034 (1.4%) 


1.834 (0.7%) 


Model III 


1.9 


1.7 


0.6802 (0.7%) 


3.466 (0.9%) 


1.083 (1.0%) 



but is related to pp and hence position-dependent. Unless stated otherwise, we fix the central (i.e . max imum) value of the 
field to be Hci{0) = 10^^ G, to fit with the estimate from iGlampedakis. Andersson fc SamuelssonI (201l|). This corresponds 
to a volume-averaged value of Hd « 3 x 10^^ G. 



6 RESULTS 

Our results are divided into four subsections. We begin by testing the convergence properties of the code and comparing our 
results with those of previous studies. We then present results for two-fluid stars in normal MHD, focussing on the effect of 
stratification. In the third subsection we study configurations where the protons form a type-II superconductor, and in the 
fourth subsection we give approximate relations for the magnetic ellipticity of stars as a function of field strength. 



6.1 Tests of the code 

Before we present results for our NS models, we test how accurate our code is by using the virial test (|56|l to determine the 
relative error in calculating an equilibrium configuration for different numbers of grid points; see figure [T] We see that the 
error decreases quadratically with resolution, and hence is second-order convergent, as intended. 



Next, we compare our results with the three nonrotating and unmagnetised two-fluid models presented in lPrix fc Rieutord 

To do this we need to convert from their dimensionless quantities to ours, using the speed of fight c = 3.00 x 10 cm 



s ^ and nuclear density p„uc = 1 .66 x 10^^ g cm ^ . Som e algebra shows that our dimensionless polytropic constants (denoted 



by the index LAG) and those of lPrix fc RieutordI (j2002l ) (with index PR) are related by 

GrigPnuc \pnucj V km / \pnuc 



k: 



LAG 



Next we consider the conversion of dimensionless masses 

12000 



Pmax 
Pnuc 



Vkm/ 



Now using these, we can make a direct comparison between our results and those of iPrix fc RieutordI (|2002l ) 



(64) 

(65) 
see table [1] 
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Ptot 
Pmax 




Figure 2. Demonstrating tlie additional flexibility possible with a two- fluid stellar model. Left: an approximation to the Douchin-Haensel 
EOS for a fluid core and a solid crust, plotting the total fluid density ptot = Pn + Pp against radius r. In Douchin and Hacnsel, the 
core is fairly close to a polytrope with index Ncotb ~ 0.7 and the crust is similar to a different polytropc, having A^crust ~ 1-5, with a 
crust-core boundary at O.OSpmax and central proton fraction a;p(0) = 0.15. In our model we have a;p(0) = 0.15 and A^coro ~ A^n = 0.7, 
while A^crust ~ Afp = 1.3. This gives us a density of ~ 0.02pniax at the crust-core boundary, with a 2 km-thick 'crust' for a 10 km radius 
star. Right: the change in proton fraction Xp against radius. We fix A^n = 1 and look at the variation of proton fraction with A^p. In the 



unstratified case, A^r, 



1, the proton fraction is constant within the star. 



6.2 Two-fluid equilibrium conflgurations in normal MHD 

Our two-fluid formalism gives us a great deal of flexibility, but for clarity we will present results using just a few canonical 
values. For unstratifled models, we set A'n = A^p = 1; whilst when we refer to stratifled models we have taken A^n ~ 1 and 
A'p — 2 unless otherwise stated (see the following flgure for justiflcation of this choice). In addition, we have fixed the central 
value of the proton fraction 2;p(0) to be equal to 0.15 for all results, although we have found that our results are virtually 
independent of the value chosen. In many figures, we show stars whose magnetic fields are (probably) unphysically high — 
of the order of 10^*^ G. This is to emphasise effects which would be less obvious at lower values. We do, however, discuss the 
scaling of our results to more typical NS field strengths too. 

In this section we present results for neutron star models composed of superfiuid neutrons and normal-fiuid protons, 
subject to a magnetic field. As discussed above, this situation may apply to an interior region of magnetars, where it is 
conceivable that the second critical field Hc2 will be exceeded. Below this critical field the protons will be superconducting, 
and we cover this case in the following subsection. 

In figure[2]we show the additional flexibility possible in two-fluid models. We can model the effect of having a fluid core sur- 
rounded by a crust, by allowing the two fluids to terminate at different points, so that there is an outer region with only one of 
the fluid s pecies — in this case the pr otons. This is shown in the left-hand plot, where we attempt to flt the equation of state de- 



scribed in lDouchin fc Haensell (|2001l ') by adjusting parameters in our model accordingly. Two- fluid models also allow us to have 



stratiflcation and consequently a non-uniform proton fraction, as one would expect in real NSs jKaminker, Hae nsel fc Yakovlev 



200 



3); see the right-hand plot. Comparing this plot with flgure 1 from Glampedakis. AnderssonX'Tandeir([201ll ) showing pro 



ton fractions in more sophisticated models, we see that choosing A'^p = 2 produces a reasonable-looking proton-fraction proflle. 

Before turning to the effect of a magnetic fleld on a two-fluid star, we consider the role of rotation on it. This allows 
us to make contact with models in previous work, as well as giving us some intuition about what to expect from magnetic 
flelds. In flgure[3]we show two conflgurations rotating at their Keplerian frequency. The left-hand star is an unstratifled model 
(A'n ~ Np — 1) and both fluid surfaces coincide. The interior density contours are different, however, because the central 
proton fraction a;p(0) = 0.15; were it a;p(0) = 0.5 then each species would have the same density contours too. The right-hand 
star is a stratifled model and in this case it is the outer fluid, the protons, that determine the star's breakup velocity — th e 



neutron surface is more rounded, without the characteristic cusped-shape of a fluid at breakup (jPrix. Novak fc Comeill2005l ). 
If we had A'^n > A'^p, the neutron fluid would be outermost; however, we do not consider this case. Instead the protons in our 
models are always outermost, for rotating conflgurations but also magnetised ones, as we will see later in this section. 

Whilst our formalism is able to deal with stars that are both highly magnetised and rapidly rotating, we will concentrate 
on non-rotating models from now on. This is because we are chiefly concerned with the effect of stratification on magnetic 
fields in stars; allowing for the extra effect of rotation obfuscates the picture. This makes our results most directly applicable 
to magnetars (which rotate very slowly), but the ellipticity formulae presented at the end of this results section are equally 
applicable to the magnetic distortions in less- magnetised NSs (like pulsars). 

We begin by looking at stars with a purely poloidal magnetic field - figure U We show contours of the streamfunction u, 
which are parallel to magnetic field lines, shading the closed-field line region and marking the 'neutral line' (where the field 
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Figure 3. Density contours of two unmagnetised stellar configurations rotating at breakup velocity, left: our canonical unstratified model 
(A^n = Np = 1), right: our canonical stratified model (A^n = 1 and A^p = 2). The bold lines show the fluid surfaces — which are coincident 
for the unstratified star. For the stratified star the neutron fluid is not rotating at breakup velocity, but the proton fluid is. 




0.2 0.4 0.6 0.8 1 1.2 1.4 




0.2 0.4 0.6 



Figure 4. Contours of the streamfunction u for a star with a purely poloidal field, corresponding to magnetic field lines. The thick arc 
indicates the stellar surface, at a dimensionless radius of unity; the shaded area indicates the closed-field line region, and the small circle 



on the X-axis shows the location of the neutral line, where the field strength vanishes. Left: no stratification (Nn 
stratified model, with A^n = 1 and A^p = 3. The closed field line region is bigger in the stratified case. 



1), right: a 



vanishes) with a circle. The field con figuration of the unstra tified model is very similar to that of a single-fiuid polytrope with 
index A'^ = 1 (see, e.g., figure 3 from lLander fc Jones! (|2009l l'). Stratification has the effect of moving the neutral line inwards 
and increasing the volume of the closed-field fine region. Fixing A'n = 1, we find that the greater the value of A^p, the larger 
the closed-field line region; we have emphasised the effect by taking A'^p = 3 rather than our canonical choice of A^p — 2. 

The toroidal-field component of a mixed-field star is defined by a function f{u), which can only be non-zero within the 
closed-field line region, and by the coefficient a (see equation (O). Increasing a increases the percentage of magnetic energy 
stored in the toroidal-field component, and also the maximum value it attains relative to the poloidal component, but it 
decreases the volume of the torus containing t he toroidal-field component; see figure [5l This seems to limit the size of the 
toroidal component, as in the barotropic case jLander fc Jones! bood : ICiolfi. Ferrari fc Gualtieril I2OI0I ). The volume of the 
torus is increased somewhat for higher values of A^p, however. 

In figure [6] we show the variation of magnetic field strength within a typical mixed-field stratified configuration. The 
field reaches a local maximum in two places — the centre of the star (corresponding to the maximum value of the poloidal 
component) and near the equatorial surface (the toroidal-component maximum). The toroidal component vanishes at the 
stellar surface, but the poloidal component extends outside the star (not shown). The configuration shown has 4.7% toroidal 
energy, around the maximum value possible within our approach. The surface field strength is about 20% of the central value, 
about 50% of the maximum toroidal field strength and also about 50% of the volume-averaged field strength B. 

Next we turn to the effect of the magnetic field on the density distribution of the star. We begin with purely poloidal 
fields, in figure [T] As in the single-fluid barotropic case, these make the star oblate, but virtually all of the distortion is in 
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Figure 5. Effect of increasing a, the coefficient of the magnetic function /{«), within a stratified star. Contours of the toroidal field 
component are shown. As a is increased, from left to right, the percentage of toroidal field increases (1%,3% and 4.5%), with a corre- 
sponding increase in the ratio B^'^^/{B^^^ + B™^") — from 0.10 to 0.18, and finally 0.34. At the same time, however, the size of the 
toroidal-field region is decreased. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 6. Magnetic field strength within a stratified nonrotating NS with 4.7% toroidal energy. The toroidal field is contained within 
the small region near the equatorial surface, reaching its maximum at 2; 0.8. The poloidal field pervades the rest of the star, reaching a 
maximum at the centre. The surface field strength is about 20% of the central value, 50% of the maximum toroidal-component strength 
and 50% of the average value B. 



the proton fluid, with the neutron fluid remaining nearly spherical. The two fluid surfaces coincide at the pole, but at the 
equator the Lorentz force distorts the proton fluid, resulting in a single-fluid region composed only of protons. To make the 
effect noticeable, we have shown a highly magnetised model, with B ~ 10^^ G. The single-fluid region i s analogous to the case 
of a r otating stratified star with no magnetic field; see figure O As mentioned in figure [2] and also bv lPrix. Novak fc Comer 
( 20051 ). i;his single-fiuid region can be thought of as an approximation to a neutron-star crust. 

The unstratified (left) and stratified (right) plots of figure [7] are very similar; the only obvious difference is that when 
A'p = 2 the protons have a larger low density region than the A'n = 1 case (this is true in unmagnetised stars too). We have 
not included density contours for mixed-field stars, as these are very similar to those in figure[7l with the contour lines pushed 
outwards slightly in the region where the toroidal component is located. 

In figure [5] we show density contours of a stratified star with a purely toroidal magnetic field. Again, the neutron fiuid 
is virtually spherical and unaffected by the magnetic field, whilst the proton fiuid is highly distorted — into a prolate 
configuration this time. Note that whilst the innermost proton-density co ntours are noticeably prolate , the outer ones become 
virtually spherical; the same effect was seen in the (single-fluid) models of Ostriker fc Hartwickl ( 19681 ). who considered mixed 
fields but with a dominant toroidal component. Related to this effect, the surfaces of the two fluids are seen to be coincident 
(or very close to being so). The unstratified case is very similar, just with more evenly-spaced proton-density contours (as 
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Figure 7. Density contours for proton fluid (solid line) and neutron fluid (dashed) in non-rotating stars distorted by purely poloidal 
magnetic flelds. The bold lines represent the fluid surfaces, which coincide at the pole but not the equator. The proton fluid is seen to 
be more distorted, due to the magnetic field. The left and right plots show our canonical unstratificd and stratified models, respectively. 




Figure 8. Density contours for proton fluid (solid line) and neutron fluid (dashed) in a non-rotating stratified star with a purely toroidal 
magnetic field. The bold line represents the two (virtually coincident) fiuid surfaces, and hence the stellar surface. Again, the proton 
fiuid is seen to be more distorted, due to the magnetic field. 



for the previous figure) — it will be siiown in tlie following subsection, when we compare toroidal fields in normal MHD and 
superconductivity. 

In figure [5] we explore mixed-field configurations with the maximum possible toroidal component (within our models, at 
least). We use two different measures of the relative strength of the toroidal component: comparing the maximum values of the 
two field components using the ratio 5™?''/ (Stof^ + Spol") and looking at the the contribution of the toroidal component to the 
total magnetic energy, ftor/fmag- Fixing Nn = 1 as usual, we find that for increasing A'^p the energy percentage increases, but 
the relative maximum value of the toroidal component decreases. In all cases the toroidal component is smaller — especially 
in terms of ftor/fmag. Note that all configurations plotted are close to spherical (with an axis ratio of r-poie/req = 0.996). Very 
highly distorted stars (with consequently strong magnetic fields) can have slightly larger values of ftor/fmag. 

Finally in this section, we plot the dependence of ellipticity on magnetic field strength, figure[TOl The relation is quadratic, 
as expected from the single-fluid case, and does not appear to depend on the stratification of the star. We show distortions 
of a purely poloidal field star, and one with a relatively strong mixed field, 3% toroidal energy. The toroidal component 
decreases the oblateness of the star, as expected, but we are not able to generate mixed-field configurations where the toroidal 
component is strong enough to induce a prolate distortion. We are only able to produce prolate distortions when the magnetic 
field is purely toroidal; this case is covered in the next section. 
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Figure 9. Mixod-ficld configurations witii tiie maximum possible toroidal component, as a function of Np. All configurations have 
A^n = 1, are nonrotating and have axis ratios close to unity. The central proton fraction a;p(0) = 0.15 in all cases, but the results 
are virtually invariant of the value chosen. We see that in more stratified stars the percentage of toroidal energy can increase, but the 
maximum toroidal field strength decreases with respect to the maximum poloidal field strength. 
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Figure 10. The cllipticity of nonrotating unstratified stars as a function of magnetic field strength. Line (a) shows the dependence for 
a poloidal field, whilst (b) shows a mixed-field configuration with 3% toroidal energy. The toroidal component decreases the oblatcness 
of the star, as expected. There is no discernable difference in the corresponding figure for a stratified star; the cllipticity appears to be 
virtually independent of stratification. 



6.3 Superconducting two-fluid equilibria with toroidal fields 

In the previous section we investigated magnetic NS models with superfluid neutrons and normal protons. This situation 
could apply in (some) magnetars, if their interior magnetic fields exceed Hc2- Most NSs, however, are likely to be composed 
predominantly of superfiuid neutrons and type-II superconducting protons. In this section we present the first NS models that 
account for this, specialising to configurations with purely toroidal magnetic fields. 

In figure [TT] we consider the variation of magnetic-field magnitude within the interior of two different NS models. The 
left-hand plot shows a star with normal protons and hence subject to the usual Lorentz force of MHD. It is also unstratified. 
The field strength is seen to drop off fairly slowly before vanishing at the stellar surface. In the right-hand plot we show 
a toroidal magnetic field in a superconducting star, with a particular choice of toroidal-field function which we refer to as 
'^•^-superconductivity' — see section 3.2. This star is also stratified, and both this and the choice of toroidal-field function 
act to 'bury' the field deeper into the star. Our other choice for the toroidal-field function, '^^-superconductivity', produces 
configurations which are (in the unstratified case) indistinguishable from the left-hand plot. This is not surprising, since the 
expression for the toroidal field in this case is of the same form as for normal MHD (again, see section 3.2). 

In figure [12] we show density distributions for three NS models with toroidal magnetic fields. We consider very highly 
magnetised stars, which serve to emphasise the magnetic distortions. All three models are broadly similar: each has a spherical 
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Figure 11. The variation in field strength within stars with purely toroidal magnetic fields. The left-hand plot shows a unstratificd 
model in normal MHD; the right-hand plot shows f •^-superconductivity in a stratified star. Both plots have the same maximum field 
strength, for easier comparison. Stratification and superconductivity both have the effect of concentrating the field in a smaller, more 
central region of the star. 
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Figure 12. The distorting effect of a toroidal field. From left to right: normal MHD, f ^-superconductivity, (.''-superconductivity. The 
plots are all for very high field strengths (of the order of lO^'^ G), since the aim of the figure is to illustrate the different ways the density 
distribution is distorted in each case. We discuss more realistic values later. 

neutron-fluid distribution, coincident neutron and proton-fluid surfaces, and prolate distortions of the proton fluid. Comparing 
the density contours for each model, we see that the superconducting models have less distortion to the proton fluid in the 
outer region, and more in the centre. The central panel (("^-superconductivity) has the unusual feature of cusps in its proton- 
fluid distribution around the pole. We show unstratifled models, but the only obvious effect of stratiflcation is that the spacing 
of the proton-fluid density contours is altered, as seen earlier in figure [T] 

We mentioned that figure [12] shows models with extremely strong magnetic fields, and we will show models with similar 
field strengths in the next figure too. Whilst these high fields are useful for emphasising certain features, the main reason for 
using them comes from our numerical method, which calculates configurations based on a given distortion (axis ratio). Since 
magnetic distortions are typically very small, this means that even the axis ratio for the least non-spherical star we can specify 
(constrained by the grid resolution) corresponds to a very strong magnetic field. This is not necessarily a problem for normal 
MHD, but superconductivity will be broken at these field strengths, which exceed the second critical field of Hc2 ~ lO^*" 
G. We are able to produce these models because the destruction of superconductivity is not built into them; the equations 
may be solved for any field strength. Having done so, however, we need to check that these models are consistent with our 
expectations for NSs with superconducting protons. If they are, we may extrapolate our results back to more realistic models, 
where B < Hc2. 

We perform this sanity check next, in figure [I3l We plot the scaling of ellipticity with average field strength B for NS 
models with toroidal magnetic fields, for normal (left) and type-II superconducting protons (right). Note that all ellipticities 
are negative, since these configurations are prolate. In the left-hand plot we find that the measured numerical values agree 
very well with the expected scaling e oc B^, with small deviations when B > 10^^ G. 

We now turn to the right-hand plot of figure [13] This time we plot the ellipticity against B, not B^ . First we look at 
configurations with (^^-superconductivity, for a central critical field value Hci{0) — 10^^ G (the points on the line marked 
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Figure 13. Left: in normal MHD the cUipticity scales with B^; right: in superconductivity it scales with H^iB - this is illustrated by 
showing the results for a central critical-field value -ffcl(O) = IQi® G (a), and also -ffci(O) = 2 X IQl^ G (c); the latter line has double the 
gradient, (a) and (c) are for (^^-superconductivity; (b) is for (^^-superconductivity with i?cl(0) = lOi® G. 



(a)) and also for Hci{Q) = 2 x 10^'' G (the points on line (c)). In both cases we see that the points lie virtually on straight 
lines, showing that e ex -B. In addition, line (c) has twice the gradient of line (a). This confirms that despite the high field 
strengths we are obliged to use (see the discussion above), we find the correct scaling of the ellipticity: e ex Hc\B. This gives 
us more confidence about our results and means we can safely extrapolate to more typical NS field strengths. We also present 
ellipticities for (^^-superconductivity (points along line (b)). The level of distortion in this case is very similar to that in the 
C,'^ case. 

Note that for both plots in figure [T3] we consider unstratified models, but we find that stratification makes no discernable 
difference to the ellipticity results. This is the same as for the previous plot showing ellipticities, figure [TOl 



6.4 Ellipticity formulae 

In figures [10] and [13] we plotted ellipticities of two-fluid stars as a function of fleld strength for poloidal, toroidal and mixed 
fields in normal MHD, and toroidal fields in superconductivity. In all cases we found the scaling of points on the plots was in 
convincing agreement with the expected results: e oc in the normal case and e cx HciB in superconductivity. In addition, 
for the very highest field strengths we were able to see slight deviations from this linear result. From the lines fitted to our 
data points we get quantitative relations between ellipticity and magnetic field strength, which we present here scaled to 
somewhat lower field strengths. In all cases the formulae are for two-fiuid models without stratification, but we found that 
stratified stars obeyed the same relations, to a good approximation. 

We begin with results for stars with superfiuid neutrons and normal protons. In the case of purely toroidal fields we have 

(66) 



1015 G 



^ = 4.3xl0-(-^) . (67) 



whilst for purely poloidal fields the relation is 

1015 G, 

We also consider mixed fields with 3% toroidal energy. This is relatively strong for our unstratified models, but still only 
produces a 14% reduction in oblateness from the purely poloidal case. 

Next we consider models comprising superfluid neutrons and type-II superconducting protons with purely toroidal mag- 
netic fields. As discussed in section 3.1, we have some flexibility in choosing the function rj which governs the magnetic fleld. 
For our flrst choice of function, (^^-superconductivity, we have 

whilst for (;^^-superconductivity the relation is 



1016 G / V 1015 G 



We recall that these results were obtained using a density-dependent first critical field Hci = hcpp, where he is a constant; 
in the above expressions we use the central critical fleld value _ffci(0), normalised to 10^*^ G. The equivalent volume- averaged 
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Figure 14. Below the (volume-averaged) second critical field Hc2 a NS has typc-II superconducting protons and its ellipticity scales 
with HciB — considerably higher than the normal-MHD curve plotted below it. Above the critical field superconductivity is destroyed 
and the ellipticity has the normal-MHD scaling, quadratic in B. 



value is Hci « 3 x 10^^ G. We note that the functional form of rj has little effect on the ellipticity relat i on. T he results in 
the superconducting case are in very good agreement with the barotropic study of Akgiin fc WassermanI ( 20081 ): using their 
ellipticity formula (71) gives a result just 5% different from our ^^-superconductivity relation (|68p and 15% different from 
l|69l). 

We summarise these last two extrapolated formulae in figure [T4l where we show the expected variation of ellipticity with 
field strength in a two-fluid NS model with a toroidal magnetic field. We assume that below the volume-averaged value of 
the second critical field Hc2 ~ 10^^ G the protons are superconducting, whilst above it they are normal. The neutrons are a 
superfiuid in both cases. In the superconducting regime, e oc B and the star is more distorted than one would expect from 
normal-MHD models. As the field strength is increased superconductivity is broken and the protons obey normal MHD, with 
eocB\ 

Note that superconductivity can make a huge difference to the ellipticities predicted for most pulsars ( Wasserman 20031 ) : 
the discrepancy from normal MHD scales with Hci/B. Let us take a pulsar with a volume-averag ed field B = 10^^ G, and 
assume this field is purely toroidal (since our only superconducting models are for this case). If the protons form a normal 
fiuid we may apply equation (|66[l to predict that its ellipticity will be —3.0 x 10"^"; if they are a type-II superconductor then 
equation (|68[) predicts that the ellipticity of the pulsar will be —2.3 x 10~^ — a factor of around 800 larger. 



7 DISCUSSION 



This paper, together with a companion paper iGlampedakis. Andersson fc Landed (|201ll ). presents the first results on the 
equilibria of multifluid neutron stars with magnetic fields. Although it has been thought for decades that the fiuid interiors of 
all but the youngest neutron stars are likely to comprise (predominantly) a neutron superfiuid and a magnetised proton fluid, 
no previous studies have constructed equilibria of this kind. A related problem is the role of stratification in NSs. This renders 
the stella r matter non-baro tropic and prevents the construction of single-fluid equilibria using the usual Grad-Shafranov 
equation ( Reisenegger 20091 ). Virtually no studies have attempted to confront this extra difficulty; to our knowledge the only 
previous studies of magnetic e quilibria in stratified stars are by Braithwaite — although these adopt an ide al-gas EOS more 
suited to main-sequence stars ( Braithwaite fc Nordlundl 20061 : Braithwaite 20091 ) — and the recent paper by Mastrano at al. 



( 201ll ). i^inally, in most neutron stars the protons are likely to form a type-II superconductor, which significantly changes the 
form of the magnetic force from the familiar Lorentz forc e. Once again, this effe c t has been neglected in most studies; the only 
related work on superconducting NS equilibria is that of I Akgiin fc Wasserman (|2008l ). although these are single-fluid models 
which do not have a separate neutron superfluid. 

Within this paper we explore some of these issues. In the normal-MHD case, we show how the equations for a magnetised 
barotropic star may be generalised quite straightforwardly to a two-fluid model in which each fluid species obeys its own 
polytropic relation. By choosing different proton and neutron polytropic indices, we have a natural way of introducing 
stratiflcation into the models. For simplicity we assume that all the stellar matter is multifluid, whereas NSs will also have 
regions which behave as a single fluid; when purely multifluid stars are better understood, it will be a logical future step to 
account for the effect of having these different fluid regions. 
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The equilibrium equations governing our stellar equilibria feature a number of apparently arbitrary magnetic functions, 
related to the strength of the magnetic field and its poloidal and toroidal field components. We argue that there are actually 
a number of restrictions on these functions, with little freedom in choosing them — and hence limited scope for producing 
qualitatively different equilibria from those we present here. Given this, we believe our results are quite generic to multifluid 
stars. 

We also discus s the equations governing equilibria wi t h typ e-II superconducting protons, starting from the relations 
derived in detail by Glampedakis. Andersson fc SamuelssonI (2011). The main difference from the normal-MHD two- fluid case 
is in the form of the magnetic force, which depends on the usual magnetic field B but also on the first critical field Hd, 
related to the presence of quantis ed fluxtubes. We use an approach inspired by the derivation of the Grad-Shafranov equation 
(see, e.g.; 



tne presence oi quantis ed nuxtubes. We use an approacn mspired by tne derivation oi tne t^rad-bnatranov equation 
Lander fc JoiieJ (|2009l )) to try to produce equilibrium equations that may be readily solved. As in the normal- 



MHD case, we arrive at a general equation governing the magnetic field, at which point one must specialise to different field 
configurations to get a solution. We remark that one important step from the normal-MHD derivation does not carry through 
to the superconducting case (and so may lead to significant differences between the two states), but we postpone a detailed 
discussion of this point to a later paper. For a purely toroidal magnetic field however, a solution may be obtained quite easily, 
so we specialise to this case in this paper. 

In all cases, the eq uilibrium equa tions are solved numerically using an iterative procedure, which is a variation of the 
Hachisu SCF method (|Hachisul llQsi ) . This metho d is quite versatile and allows us t o con sider extremely strong mag- 
netic fields, in contrast with the study reported in iGlampedakis. Andersson fc Landed (|201ll ). where the magnetic field is 
regarded as a perturbation on a spheric al background star. The good agreement between the results of this paper and 



Glampedakis. Andersson fc Landed ((20111), however, vindicate the latter work's perturbative ansatz. 



One main aim of this paper is to explore the role of stratification in a magnetised neutron star. We increase the stratifi- 
cation of the star by increasing the value of the proton-fluid polytropic index A^p, whilst fixing A^'n = 1. For a purely poloidal 
field, this results in a larger closed-field line region, with the 'neutral line' (where the field strength vanishes) moving inwards; 
see figure m In mixed-field configurations the toroidal field component is contained in this region, and is consequently larger 
in stratified stars. Even in the stratified case, however, the energy in the poloidal component is always considerably greater 
than that in the toroidal component. This is because the parameter which increases the toroidal-component strength also 
decreases the size of the region occupied by this component (figure [5]). These results are in good agreement with those of 



Glampedakis. Andersson fc Landed (|201ll ). Finally, for magnetic ellipticity relations our results suggest that stratification is 



not very important; there is little difference from the single-fiuid barotropic relations. 

We have not investigated the stabili ty of any configurations p resented here, although purely toroidal fields are unstable 
in both stratified and unstratified stars ( Goossens fc Tavled ngSOf), and it seems likely that the same will be true for purely 
poloidal fields, as in the unstratified case (|Markev fc Tavled Il973l l. This leaves mixed poloidal-toroidal fields as the only 
candidates for stable configurations, and hence the best models for NS magnetic fields. Although our configurations have 
quite small percentages of magnetic energy in the toroidal component, due to the small volume of the star it occup ies, the 
magn itude of the two field compone nts can be comparable. This suggests that these models could indeed be stable (jWright 
1973h . Note that ISraithwaitel (|2009h suggested that a considerably higher percentage of toroidal-field energy was needed for 
stability. That work considered very different stellar models, however, more applicable to main-sequence stars than neutron 
stars. There is no consensus, as yet, about what percentages of each field component are likely to exist in real NSs. 

Our other results concern equilibrium configurations of neutron stars with superfluid neutrons and type-II superconducting 
protons. In thi s initial study we specialise to purely toroidal magnetic fields, even though the equilibria we generate are likely 
to be unstable ( Akeciin fc Wasserman|[200^ . From our results we are able to produce the first formulae for magnetic ellipticities 
of superconducti ng NSs based on a two-fluid model. Our results suggest that the distortion is around 50% greater than the 
estimate used bv lCutled (|2002l ) (having accounted for the fact that we also use a larger value for the critical field Hci)- There 
is very good agreement, however, with the ellipticity formula (equation 71) from Akgiin fc Wasserman ( 20081 ) . If a substantial 
proportion of prot ons in a typical NS inter ior form a type- I I sup erconductor, then ellipticity formulae based on models in 
normal MHD (e.g. Haskell et al. ( 20081 ') and Lander fc Jones ( 20091 )') greatly underestimate the potential magnetic distortion. 

With the two-fiuid model discussed in this paper, we are able to produce a wide range of different NS equilibria, accounting 
for stratification, superconductivity (with toroidal fields) and rotation. As well as being interesting in their own right, these 
equilibria could be used as background models on which to study perturbations. This would help us understand the stability 
of these various models, and their oscillation modes. The next logical step, however, is to consider superconducting stars with 
purely poloidal and mixed poloidal-toroidal magnetic fields. We hope to tackle this problem in a future study. 
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